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Abstract 

We perform a non-perturbative determination of the improvement coefficient ca 
to remove 0{a) discretization errors in the axial vector current in three- flavor lattice 
QCD with the Iwasaki gauge action and the standard 0(a)-improved Wilson quark 
action. An improvement condition with a good sensitivity to ca is imposed at constant 
physics. Combining our results with the perturbative expansion, ca is now known 
rather precisely for a~ l >l.Q GeV. 
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1 Introduction 



Various discretizations of QCD on a lattice are presently used in the large scale efforts 
aiming at non-perturbative results in the theory of strong interactions (see [U El El HI El 
El El El HO] and references therein). Wilson's original formulation [TT] is theoretically very 
well founded [T2"[ [T3] and rather simple to implement in numerical simulations. The flavor 
symmetries are exact and with modern algorithms [HI [T5], [16j [TTJ, [18] the regime of small 
quark masses and small lattice spacings can be reached [H [19]. On the other hand it is 
well known that since the chiral symmetries are broken by the Wilson term, lattice artifacts 
linear in the lattice spacing are present. It has long been understood how these can be 
removed by applying Symanzik's improvement programme [201 EH E21 E31 El]- One has to 
add dimension five fields to the lattice Lagrangian and (for example) dimension four fields 
to the quark bilinears. In particular the bare flavor axial current 

A;(x)=^(x)T a ^ l5 ^(x) (1) 

(the SU(iVf) generator T a acts in flavor space) is improved by 

(A)*(x) = Al{x)+ac K l -{d, + d;)P a {x), P a (x)=^(x)T a 75 ij(x) } (2) 
with , 1 

d,f(x) = -[/(* + a/2) -/(*)], % f(x) = ~[f(x) - f(x - aft)] . (3) 

The coefficients of these correction terms, such as ca, can be determined non-perturbatively 
by requiring specific continuum chiral Ward-Takahashi identities to be valid at finite lattice 
spacing [25J. One is then sure that the 0(a) effects are entirely removed. Details of this 
programme as well as the present status have recently been reviewed [26J. Here we just 
mention that the coefficient c= w of the Sheikholeslami- Wohlert term [21] , the only dimension 
five correction to the actional, has been determined non-perturbatively for different gauge 
actions and numbers of flavors [251 EZ1 EH1 EH] ■ 

Next to c sw , the axial current improvement coefficient ca is of particular relevance - for 
example in the determination of weak leptonic decay constants such as F n or the quark 
masses. Non-perturbative determinations of ca have been studied for Nf — and 2 in 
refs. [251 EH Ell E21 E3] ■ It turned out that they need special care since the spread between 
CA-values computed from different improvement conditions is significant around a « 0.1 fm. 
There is nothing fundamentally wrong with this fact. However, as explained in some detail 
in refs. [2E1 Ell EH El] , in such a situation it is important to impose improvement conditions 
on a line of constant physics. This means that as the lattice spacing a is varied, all other 
physical scales are kept fixed. The remaining effects (after improvement) are then smooth 
0(a 2 ) terms. 

Here we apply this strategy to the theory with iVf = 3 flavors and the Iwasaki gauge action 
[33] , which is of immediate interest to the large scale computations of the CP-PACS/JLQCD 

1 We neglect small 0(am) modifications of the gauge couplings and quark masses [52 |23], as they are not 
so relevant in practice [52] . Here m stands for any of the quark masses. 
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Collaborations [3] . All known practical methods for a computation of c\ start from the fact 
that in the continuum limit the (PCAC) quark mass 

m - 2(0'|P«|0) [ } 

does not depend on the choice of the external states \<f>) , \4>'). This is just a rephrasing 
of the PCAC (operator) identity. On the lattice an 0(a) dependence will exist in general. 
It is reduced to 0(a 2 ) by improvement. Requiring m to be the same for two different 
choices of \<p) , | </>'), or as we will say later "two different kinematical conditions", allows a 
determination of ca when c sw is already known. 

As in ref. [33], we use the Schrodinger functional defined in a Euclidean L 3 x T world to 
construct suitable states with a large sensitivity to ca- In the following section we define the 
exact choices of kinematical conditions. The reader who is familiar with ref. [33] may skip 
this section and proceed directly to the description of the simulation details, sect. El and the 
results, sect. HI We finish with some conclusions. 



2 Improvement condition 

We introduce the following Schrodinger functional [HHl EZ] correlation functions 



ft\x ) = -yJ>S(*)0 a >M>, f^(xo) = -jJ2( pa ^° a ' in) ) and ( 5 ) 

X X 

j(n,m) _ _±(^Qla,(n)Qa,{m)^ ^ 



with 



a,(n) = a £y») (y _ z) f (y) T a 75 C(z). (7) 



where ( and ( are the fermionic boundary fields on the Xq = timeslice (O' is defined in 
the same way in terms of the boundary fields at x = T). The correlators depend on the 
smooth functions cu^. Here, as in ref. [33], we use three wave functions 



,(») = 

kez 3 



cu (1) (r) oc e-^ ao , c^(r) oc (|r|/r ) e-^ aa , cu (3) (r) oc e~^^ , (9) 

with ao — L/6. The normalization factors jV( n ) are fixed by a 3 ^ x (lo>^) 2 = 1. 

By suitably combining the operators O a ^ n \ the resulting correlation functions get contri- 
butions from different states in the pseudoscalar channel. In fact we construct the boundary 
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operators O and Ox, which mainly couple to the ground and first excited states respectively, 
by using the eigenvectors of the 3x3 symmetric matrix f^ l,m ^ 

O- = J2r] ( n) O a ' {n \ O a 1 = ^r 1 ^ ) 0< n \ (10) 

n n 

where r] (771) represents the eigenvector associated with the largest (2nd largest) eigenvalue. 
The corresponding correlators fx,i = ^2 n Vi fx with £ = 0,1 and X=A,P are eventually 
used to define the improvement condition, which reads 

m (x ) =m 1 (x ) 1 (11) 

where 

mj(xo) = r i (x ) + c A as i (x ), (12) 
Solving eq. (TTTT) for ca yields 

1 Af ( Tn ) 

ca(^o) = — "t ~i — c ) Ar(x ) = ri(x ) - r (x ), As(x ) = Si(ar ) - s (x ) . (14) 
a As{x ) 

The sensitivity of the improvement condition to ca is given by a |As(x )|- In the ideal case of 
exact projection on the ground (ir) and first excited (n 1 ) states (and large T — x ) that would 
be given by a {m\ — As discussed in ref. [33] however the vectors 77$ do not achieve 
perfect projection and the correlator /a, 1(^0) for example gets some contribution from the 
ground state. Anyway what is really needed is that at intermediate times xq — T/2 where 
we extract ca, the correlation functions are dominated by states with different energies, such 
that the sensitivity is high. We will see in sect. H]that in our setup this is indeed the case. 



3 Simulation details 

We work in the theory with three (dynamical) degenerate flavors of non-perturbatively im- 
proved Wilson-fermions and the Iwasaki gauge action [2U, EHJ [35] . The latter reads 

S 9 = A Yl <>o)^ReTr[l-P^(a;)]+ ^ w* u (x )±ReTr[l-R^(x)]\ , (15) 

I. X,[l<V X,fJ,,U J 

where (3 = and P^ u and R^ u are the lxl and 1x2 Wilson loops in the (/i, v) plane. 
Their weights are w^ u = 3.648 and w^ u = — 0.331 on periodic lattices. 



4 



The Schrodinger functional formalism is implemented on a L 3 x T lattice with L = T. 
The background field is set to zero and the fields are chosen to be periodic in space. The 
weights in the gauge action are modified to the following choice [38J 

p, v _ J 1/2 at x = or T, and /i, u^4 , . 

w ^ x °> = \ 3.648 otherwise { } 

{0 at x = or T, and fi, v ^ 4 

-0.331 x (3/2) at x = or T, and /i = 4 (17) 
—0.331 otherwise 

which entails tree level 0(a) improvement "at the boundaries" [22]- The coefficients of the 
O(a) boundary counterterms for the fermions are also set to their tree level value [22] • Note 
that this is not at all essential. Irrespective of whether the boundary improvement terms 
are implemented, eq. (lllj) is a correct improvement condition |22j . 

We simulate at three points in the ((3, L/a, k) space on a line of constant physics defined 
by keeping the volume and the quark mass fixed. Scales are fixed through r [39j and we will 
use r = 0.5 fm to quote physical units. For our action, the ratio r /a has been computed 
in the region 1.83 < (3 < 2.05 [3"ffi With a slight interpolation of the data of ref. [3] we 
fixed L/r « 3, somewhat larger than the physical size used in ref. [33]. The resulting pairs 
(P,L/a), together with some algorithmic details are collected in table [H 



p 


L/a 


K 


^MD 


-Wpoly 


iVtraj 


-PpHMC 


-Pnmt 


1.83 


12 


0.13852 


90 


200 


3800 


0.90 


0.97 


1.83 


12 


0.13867 


90 


220 


3800 


0.89 


0.95 


1.95 


16 


0.13685 


125 


230 


3000 


0.91 


0.96 


1.95 


16 


0.13697 


140 


260 


3000 


0.94 


0.94 


2.05 


20 


0.13604 


130 


350 


3000 


0.87 


0.97 



Table 1: Simulation parameters. We denote the number of the Molecular Dynamics steps by Nmd, 
the order of the polynomial approximation by N po \ y and the number of unit length trajectories by 
-^traj • The acceptance rates for the PHMC updating and for the noisy Metropolis test in the PHMC 
algorithm are denoted by -Pphmc a nd -Pnmt respectively. 



The hopping parameter k is tuned in order to give a bare quark mass m re f of about 
15 MeV. At /3= 1.83 and 1.95, two quark masses around 10 and 20 MeV are simulated so 
that we can interpolate ca to m re f. Notice that we are ignoring presumably small changes 
of the renormalization factors in our range of (3 and we just keep the bare quark mass fixed. 
Also, the 1-loop value of c\ [38], HQ] is used at this point in the definition of the quark mass. 

The algorithm has been described in ref. [TS]. It is a combination of HMC [41J and 
PHMC H3J . Non-Hermitian Chebyshev polynomials P[D] are used to approximate the 
inverse square root of the Dirac operator D required for the third flavor, whereas the other 

2 In ref. [H] r$/a is extrapolated to the physical point, with the strange quark mass determined from the 
physical mass of the Kaon. 
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two flavors are treated using the usual HMC pseudo-fermion action. The number of molecular 
dynamics steps is chosen such that the acceptance rate Pphmc is about 90%. In order to 
make this algorithm exact, the correction factor 



is taken into account by a noisy Metropolis test [Hj. The order of the polynomial approxi- 
mation is chosen such that the acceptance rate of the noisy Metropolis test is around 95%. 
Throughout all the computation the symmetrically even-odd preconditioned version of the 
Dirac operator [15], H5] is used. 

The correlators in eqs. (0, ED are measured each 5th trajectory and residual autocorrela- 
tions are estimated by binning the jack-knife samples. For ca(xq) the errors flatten out for 
bin-sizes larger than four, which is what we finally use in our analysis. 

4 Numerical results 

4.1 Wave function projection 

As discussed above, the analysis starts with the determination of the eigenvectors of the 
correlator matrix /-["'"^(n, m = 1,2,3). The results at the lightest quark mass for each (3 
value are given in table [21 The errors on the components of the eigenvectors are less than 
10 -3 and 10~ 2 for i]q and rji, respectively. 

The entries of the (normalized) eigenvectors are ratios of correlation functions for which 
the Z-factors of the boundary fields cancel. They will thus have a well defined continuum 
limit along a line of constant physics as long as the wave functions a/ n ) are defined only 
through physical length scales. Indeed, we observe only a small lattice spacing dependence 
and also note that the values we obtain are close to those in ref. |33j, where slightly smaller 
values of do (~0.2 fm) and L are used in two-flavor QCD with the plaquette gauge action. 

Using the measured eigenvectors, we now construct the pseudoscalar correlators /x,o an d 
/x,i (X = P, A). Fig. [T] shows the effective masses in units of the box size for the projected 
correlators /p^ at /3 = 1.83 and 1.95. Clearly, the correlators are dominated by different 
states and the effective masses are well separated even for large times. The data at the two 
coarser lattice spacings are obtained at physical quark masses similar to each other and we 



P 



corr 



det [W[D]] , W[D] = P[D]D 



(18) 



(3 K 




1.83 0.13867 
1.95 0.13697 
2.05 0.13604 



( 0.5459, 0.5920, 0.5929 ) ( 0.8323, -0.3019, -0.4649 ) 
( 0.5415, 0.5942, 0.5948 ) ( 0.8367, -0.312 , -0.4500 ) 
( 0.5360, 0.5962, 0.5976 ) ( 0.8371, -0.2836, -0.4679 ) 



Table 2: Example of eigenvectors r]Q and r\\ in three-flavor QCD at each f3. 
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observe good agreement also for the effective masses in the pseudoscalar channel (mp L m 3 
and m-p iL m 11). In table [3] we have included also the combination L 2 As, which has a 
continuum limit when all quantities are computed on a line of constant physics. While an 
a-dependence appears to be present in this combination, this is small. We can take its 
smallness as good evidence that our improvement condition does not suffer from large 0(a 2 ) 
contributions. 



20 



-J 



15 



10 



5- 



P= 1.83, K = 0.13867 



e e 



e © 



o / PJ 
□ f P 



e © 



20 ~ P = 1.95, K = 0.13697 



15 



10 



5 - 



o 4, 
□ /p 



n m 



10 



e © 

_j i 

5 



e e © 



_i 

10 



15 



Fi gure 1: Effective mass mp e g — ^ ^°s{fp,i{ x o ~ a ) / fp ,i{ x o a )) f° r /p,o and fp ; i rescaled by L. 
The dotted line shows where the effective mass is equal to a . 



Note, however, that the effective mass mp^ at our smallest (3 is already close to the 
cutoff (i.e. mp ; i L ~ L/a in fig. [T]). This implies that a rapid increase of the residual 0(a 2 ) 
effects might occur if one were to evaluate the improvement condition at even coarser lattice 
spacings. 



4.2 The improvement coefficient 

With the projected correlation functions we can proceed to the extraction of c\ itself. Tabled 
lists the differences Ar and As for the lightest quark mass at each (3. In all cases we see a 
good signal for As and thus have a large sensitivity to ca- 



p 


K 


a Ar 


a 2 As 


L 2 As 


1.83 


0.13867 


0.0229(14) 


0.429(22) 


62(3) 


1.95 


0.13697 


0.0072( 7) 


0.236(14) 


60(4) 


2.05 


0.13604 


0.0036( 3) 


0.133( 6) 


53(2) 



Table 3: Examples of Ar and As at xq = T/2. 



In fig. El we plot the effective ca(xq), cf. eq. ( JT4l) . from the finest and coarsest lattices. 
The value of ca(xq) stabilizes after only a few lattice spacings from the lower temporal 
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boundary where higher exited states contribute. In all cases Xq = T/2 is already in the 
region, where these effects are small and we use this choice to complete the definition of ca- 
Results for the improvement coefficient and the PCAC quark mass from all simulations are 
collected in table HI 



0.00 r 



-0.05 



-0.10 



P=1.83 

O K= 0.13852 
□ K= 0.13867 



x Q la 



0.00 



-0.05 



-0.10 



P=2.05 

O K = 0.13604 



10 



10 

x Q /a 



o O 



15 



Figure 2: Effective value of ca as a function of xq at (3 = 1.83 (left panel) and (3 = 2.05 (right 
panel). The data at /3= 1.83 and k = 0. 13867 are slightly shifted along the vertical axis for clarity. 



p 


K 


am 


ca 


2.05 


0.13604 


0.00554(14) 


-0.0272(18) 


1.95 


0.13685 


0.01020(29) 


-0.0348(25) 




[interp.] 


am ref 


-0.0319(18) 




0.13697 


0.00508(28) 


-0.0303(24) 


1.83 


0.13852 


0.01406(54) 


-0.0519(23) 




[interp.] 


am ref 


-0.0528(17) 




0.13867 


0.00614(63) 


-0.0534(24) 



Table 4: Numerical results for am and ca- Also shown are the results of the interpolation to m re f 
at the two coarser lattice spacings. 



4.3 Interpolation of ca 

As discussed above, we aim at evaluating the improvement condition on a line of constant 
physics in order to avoid potentially large O(a) ambiguities in ca itself. To this end we 
interpolate the results for ca at (3 = 1.83 and 1.95 to a quark mass m re f that is matched to 
the one measured on the finest lattice. The quark mass dependence seems to be very small 
and thus the uncertainties in the quark masses become unimportant and we obtain ca at 
m re f with a small statistical error. 
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For future use we summarize the present results for the improvement coefficient in an 
interpolating formula (iT9"l) . which by construction reduces to the one- loop result from refs. [3E1 
HO] in the perturbative limit 



/ 9n •> 1 - 0.195 gl 

c A(9 ?) = -0.0038,? x r -^| 



(19) 



It is plotted in fig [3j where one can verify that this formula reproduces the data well and 
gives a smooth interpolation in the range of (3 values we simulated. As was the case with 
the plaquette gauge action and two quark flavors [33], the non-perturbative result is quite 
different from the one-loop estimate for practically relevant lattice spacings. 



0.00 



-0.02 



-0.04 



-0.06 - 




Figure 3: Non-perturbative estimate of ca as a function of g$. 



4.4 Systematic uncertainties 

The computation of ca on a line of constant physics reduces the intrinsic ambiguity on the 
improvement coefficient to a smooth O(a) form. Deviations from this condition will lead 
to systematic effects and we should therefore check the consequences of variations of the 
physical volume and quark mass for our improvement condition. 

All simulations, on which we report here, are performed at fixed physical volume and 
we thus have no direct check of the volume effects on ca from this improvement condition. 
From [33] we know that those can be large, but we know that our condition guarantees that 
they disappear smoothly as we approach the continuum limit, especially since our volume 
scaling is based on actual measurements of r /a. 

From the data at the two coarser lattice spacings in table H] it is evident that the quark 
mass dependence of ca is very weak in our setup. This implies that no fine tuning of 
m is required and also a posteriori justifies the fact that we ignore small changes of the 
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renormalization factor in our range of (3 and use the bare quark mass in our definition of a 
line of constant physics. 

As mentioned above, the energy of the first excited state at our lowest /? = 1.83 is close to 
a -1 . Consequently, enforcing the present improvement condition at (5 < 1.83 may induce large 
0(a 2 ) scaling violations in the axial current. While larger volumes might help in lowering 
this energy, the observation shows that even with improved gauge actions, one should not 
push the simulations too much towards coarse lattice spacings. On the other hand it is useful 
to repeat our earlier observation: within the range of lattice spacings covered here, we see a 
reasonable scaling of L 2 As; this is a good hint that the considered matrix elements do not 
suffer from large a-effects. In retrospect the same statement can be made about the N? = 2 
computation with plaquette gauge action [33J. 

5 Conclusions 

We have computed the 0(a)-improvement coefficient ca(<7o) of the axial current non-per- 
turbatively in three-flavor QCD with the Iwasaki gauge action and non-perturbative c sw (go) 
[29] . The improvement coefficient ca is parametrized as a function of go- Since the results 
connect smoothly to the one-loop formula at weak coupling, a simple interpolation formula 
(fT§|) could be given in the range of a -1 > 1.6 GeV. 

We note that at the largest lattice spacing covered, the correction term amounts to 10- 
15% in decay constants F n , Fx and then also in the renormalized quark masses evaluated 
from the PCAC relation. As a next step, a full non-perturbative evaluation of these quan- 
tities now requires the computation of the renormalization factor and for interesting 
applications of vector form factors also the corresponding quantities Zy and cy are very 
relevant. On the other hand, improvement terms proportional to the light quark masses are 
suppressed by the smallness of am. It then appears sufficient to approximate the associated 
coefficient by one- loop perturbation theory [381 00] . 
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